library(contrast)
library(data.table)
library(interflex)
library(lme4)
library(lubridate)
library(merTools)
library(stringr)
library(texreg)

# SET WORKING DIRECTORY TO LOCATION OF REPLICATION ARCHIVE
setwd("")

responses <- fread("chevron_data.csv", header = TRUE, stringsAsFactors = FALSE)

################################################################################

# TABLE 1

# H1: treated bureaucrats are less likely to stay in their jobs
likelihood_job <- lm(likelihood_worknextyear ~ chevron_condition, data = responses[which( responses$state!="Indiana" & responses$state!="Nebraska" & responses$state!="North Carolina"),])
summary(likelihood_job)

# H2: treated bureaucrats are less likely to stay in their jobs as PSM increases
likelihood_job_int <- lm(likelihood_worknextyear ~ chevron_condition*psm, data = responses[which( responses$state!="Indiana" & responses$state!="Nebraska" & responses$state!="North Carolina"),])
summary(likelihood_job_int)

# H3: treated bureaucrats are less likely to invest effort
invest_effort <- lm(invest_effort ~ chevron_condition, data = responses[which( responses$state!="Indiana" & responses$state!="Nebraska" & responses$state!="North Carolina"),])
summary(invest_effort)

# H4: treated bureaucrats are less likely to invest effort as PSM increases
invest_effort_int <- lm(invest_effort ~ chevron_condition*psm, data = responses[which( responses$state!="Indiana" & responses$state!="Nebraska" & responses$state!="North Carolina"),])
summary(invest_effort_int)

wordreg(l=list(likelihood_job, likelihood_job_int,
                invest_effort, invest_effort_int),
         file = "Chevron_mainresults.doc",
         custom.coef.map = list("(Intercept)" = "Intercept",
                                "chevron_conditionTreatment" = "Chevron Treatment",
                                "psm" = "Public Service Motivation (1-5)",
                                "chevron_conditionTreatment:psm" = "Chevron Treatment:Public Service Motivation"),
         custom.model.names = c("Likelihood of Remaining in Job", "Likelihood of Remaining in Job", 
                                "Effort Invested in Expertise", "Effort Invested in Expertise"),
         caption = "Effect of Chevron Prime on Turnover Intention and Effort Invested in Expertise",
         caption.above = TRUE,
         stars = 0.10,
         include.rsquared = FALSE,
         include.adjrs = FALSE)

################################################################################

# TABLE SI.1

# Numbers of respondents in sampling frame from each state as reported in the
# table, response rates calculated as follows:

# Connecticut
dim(responses[which(responses$state=="Connecticut")])[1]/11929

# Illinois
dim(responses[which(responses$state=="Illinois")])[1]/1840

# Indiana
dim(responses[which(responses$state=="Indiana")])[1]/27168

# Nebraska
dim(responses[which(responses$state=="Nebraska")])[1]/12701

# New Hampshire
dim(responses[which(responses$state=="New Hampshire")])[1]/8898

# North Carolina
dim(responses[which(responses$state=="North Carolina")])[1]/48283

# Oregon
dim(responses[which(responses$state=="Oregon")])[1]/36460

# Vermont
dim(responses[which(responses$state=="Vermont")])[1]/3645

# Overall
dim(responses)[1]/150924

################################################################################

# TABLE SI.2

# gender
table(responses$gender[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))])
round(prop.table(table(responses$gender[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))]))*100,1)

# age
table(responses$age[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))])
round(prop.table(table(responses$age[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))]))*100,1)

# income
table(responses$income[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))])
round(prop.table(table(responses$income[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))]))*100,1)

# education
table(responses$education[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))])
round(prop.table(table(responses$education[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))]))*100,1)

# race
table(responses$race[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))])
round(prop.table(table(responses$race[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))]))*100,1)

# Hispanic
table(responses$hispanic[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))])
round(prop.table(table(responses$hispanic[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))]))*100,1)

# PID
table(responses$pid7[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))], useNA = "always")
round(prop.table(table(responses$pid7[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))], useNA = "always"))*100,1)

# ideology
table(responses$ideology[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))])
round(prop.table(table(responses$ideology[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))]))*100,1)

# years exp total
table(responses$years_exp_total[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))])
round(prop.table(table(responses$years_exp_total[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))]))*100,1)

# appointed/civil service
table(responses$job_classification[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))])
round(prop.table(table(responses$job_classification[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))]))*100,1)

# policymaking
table(responses$jobresponsibilities_policymaking[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))])
round(prop.table(table(responses$jobresponsibilities_policymaking[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))]))*100,1)

# implementation
table(responses$jobresponsibilities_implementation[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))])
round(prop.table(table(responses$jobresponsibilities_implementation[which(!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort))]))*100,1)

################################################################################

# TABLE SI.3

# H1: treated bureaucrats are less likely to stay in their jobs

likelihood_job_all <- lm(likelihood_worknextyear ~ chevron_condition, data = responses)
summary(likelihood_job_all)

# H2: treated bureaucrats are less likely to stay in their jobs as PSM increases

likelihood_job_int_all <- lm(likelihood_worknextyear ~ chevron_condition*psm, data = responses)
summary(likelihood_job_int_all)

# H3: treated bureaucrats are less likely to invest effort

invest_effort_all <- lm(invest_effort ~ chevron_condition, data = responses)
summary(invest_effort_all)

# H4: treated bureaucrats are less likely to invest effort as PSM increases

invest_effort_int_all <- lm(invest_effort ~ chevron_condition*psm, data = responses)
summary(invest_effort_int_all)

wordreg(l=list(likelihood_job_all, likelihood_job_int_all,
               invest_effort_all, invest_effort_int_all),
        file = "Tables\\ChevronExp_allstates.doc",
        custom.coef.map = list("(Intercept)" = "Intercept",
                               "chevron_conditionTreatment" = "Chevron Treatment",
                               "psm" = "Public Service Motivation (1-5)",
                               "chevron_conditionTreatment:psm" = "Chevron Treatment:Public Service Motivation"),
        custom.model.names = c("Likelihood of Remaining in Job", "Likelihood of Remaining in Job", 
                               "Effort Invested in Expertise", "Effort Invested in Expertise"),
        caption = "Effect of Chevron Prime on Turnover Intention and Effort Invested in Expertise",
        caption.above = TRUE,
        stars = 0.10,
        include.rsquared = FALSE,
        include.adjrs = FALSE)



################################################################################

# TABLE SI.4

# CREATE NEW DATA FRAME FOR ONLY RESPONDENTS WHO ENGAGE IN POLICYMAKING SOMEWHAT
# OR VERY FREQUENTLY
responses_policymaking_often <- responses[which(responses$jobresponsibilities_policymaking=="Somewhat frequently" |
                                                  responses$jobresponsibilities_policymaking=="Very frequently"),]

# H1: treated bureaucrats are less likely to stay in their jobs
likelihood_job <- lm(likelihood_worknextyear ~ chevron_condition, data = responses_policymaking_often[which(responses_policymaking_often$state!="Indiana" & responses_policymaking_often$state!="Nebraska" & responses_policymaking_often$state!="North Carolina"),])
summary(likelihood_job)

# H2: treated bureaucrats are less likely to stay in their jobs as PSM increases
likelihood_job_int <- lm(likelihood_worknextyear ~ chevron_condition*psm, data = responses_policymaking_often[which(responses_policymaking_often$state!="Indiana" & responses_policymaking_often$state!="Nebraska" & responses_policymaking_often$state!="North Carolina"),])
summary(likelihood_job_int)

# H3: treated bureaucrats are less likely to invest effort
invest_effort <- lm(invest_effort ~ chevron_condition, data = responses_policymaking_often[which(responses_policymaking_often$state!="Indiana" & responses_policymaking_often$state!="Nebraska" & responses_policymaking_often$state!="North Carolina"),])
summary(invest_effort)

# H4: treated bureaucrats are less likely to invest effort as PSM increases
invest_effort_int <- lm(invest_effort ~ chevron_condition*psm, data = responses_policymaking_often[which(responses_policymaking_often$state!="Indiana" & responses_policymaking_often$state!="Nebraska" & responses_policymaking_often$state!="North Carolina"),])
summary(invest_effort_int)

wordreg(l=list(likelihood_job, likelihood_job_int,
               invest_effort, invest_effort_int),
        file = "Tables\\ChevronExp_constrained_policymaking.doc",
        custom.coef.map = list("(Intercept)" = "Intercept",
                               "chevron_conditionTreatment" = "Chevron Treatment",
                               "psm" = "Public Service Motivation (1-5)",
                               "chevron_conditionTreatment:psm" = "Chevron Treatment:Public Service Motivation"),
        custom.model.names = c("Likelihood of Remaining in Job", "Likelihood of Remaining in Job", 
                               "Effort Invested in Expertise", "Effort Invested in Expertise"),
        caption = "Effect of Chevron Prime on Turnover Intention and Effort Invested in Expertise",
        caption.above = TRUE,
        stars = 0.10,
        include.rsquared = FALSE,
        include.adjrs = FALSE)

################################################################################

# TABLE SI.5

# CREATE NEW DATA FRAME FOR ONLY RESPONDENTS WHO ENGAGE IN IMPLEMENTATION SOMEWHAT
# OR VERY FREQUENTLY

responses_implementation_often <- responses[which(responses$jobresponsibilities_implementation=="Somewhat frequently" |
                                                    responses$jobresponsibilities_implementation=="Very frequently"),]

# H1: treated bureaucrats are less likely to stay in their jobs

likelihood_job <- lm(likelihood_worknextyear ~ chevron_condition, data = responses_implementation_often[which(responses_implementation_often$state!="Florida" & responses_implementation_often$state!="Indiana" & responses_implementation_often$state!="Nebraska" & responses_implementation_often$state!="North Carolina"),])
summary(likelihood_job)

# H2: treated bureaucrats are less likely to stay in their jobs as PSM increases

likelihood_job_int <- lm(likelihood_worknextyear ~ chevron_condition*psm, data = responses_implementation_often[which(responses_implementation_often$state!="Florida" & responses_implementation_often$state!="Indiana" & responses_implementation_often$state!="Nebraska" & responses_implementation_often$state!="North Carolina"),])
summary(likelihood_job_int)

# H3: treated bureaucrats are less likely to invest effort

invest_effort <- lm(invest_effort ~ chevron_condition, data = responses_implementation_often[which(responses_implementation_often$state!="Florida" & responses_implementation_often$state!="Indiana" & responses_implementation_often$state!="Nebraska" & responses_implementation_often$state!="North Carolina"),])
summary(invest_effort)

# H4: treated bureaucrats are less likely to invest effort as PSM increases

invest_effort_int <- lm(invest_effort ~ chevron_condition*psm, data = responses_implementation_often[which(responses_implementation_often$state!="Florida" & responses_implementation_often$state!="Indiana" & responses_implementation_often$state!="Nebraska" & responses_implementation_often$state!="North Carolina"),])
summary(invest_effort_int)

wordreg(l=list(likelihood_job, likelihood_job_int,
               invest_effort, invest_effort_int),
        file = "Tables\\ChevronExp_implementation.doc",
        custom.coef.map = list("(Intercept)" = "Intercept",
                               "chevron_conditionTreatment" = "Chevron Treatment",
                               "psm" = "Public Service Motivation (1-5)",
                               "chevron_conditionTreatment:psm" = "Chevron Treatment:Public Service Motivation"),
        custom.model.names = c("Likelihood of Remaining in Job", "Likelihood of Remaining in Job", 
                               "Effort Invested in Expertise", "Effort Invested in Expertise"),
        caption = "Effect of Chevron Prime on Turnover Intention and Effort Invested in Expertise",
        caption.above = TRUE,
        stars = 0.10,
        include.rsquared = FALSE,
        include.adjrs = FALSE)

################################################################################

# SECTION C.3

# In-text mention of central tendency of PSM across respondents
summary(responses$psm[which((!is.na(responses$likelihood_worknextyear)|!is.na(responses$invest_effort)))])

# FIGURE SI.1
interflex(Y = "invest_effort", D = "chevron_condition", X = "psm", 
          data = responses[which( responses$state!="Indiana" & responses$state!="Nebraska" & responses$state!="North Carolina"),], 
          estimator = "raw", vcov.type = "robust", main = "Marginal Effects", 
          Xlabel = "Public Service Motivation", Ylabel = "Effort Invested in Expertise",
          ylim = c(-1, 1), na.rm = TRUE, theme.bw = TRUE, ncols = 3)

# FIGURE SI.2
interflex(Y = "invest_effort", D = "chevron_condition", X = "psm", 
          data = responses[which( responses$state!="Indiana" & responses$state!="Nebraska" & responses$state!="North Carolina"),], 
          estimator = "binning", vcov.type = "robust", main = "Binning Estimator",
          Xlabel = "Public Service Motivation", Ylabel = "Effort Invested in Expertise",
          Dlabel = "Chevron Prime",
          ylim = c(-1, 1), na.rm = TRUE, theme.bw = TRUE)

# KERNEL ESTIMATOR ANALYSIS MENTIONED IN FOOTNOTE 9 OF SI
interflex(Y = "invest_effort", D = "chevron_condition", X = "psm", data = responses[which( responses$state!="Indiana" & responses$state!="Nebraska" & responses$state!="North Carolina"),], estimator = "kernel", vcov.type = "robust", main = "Marginal Effects", ylim = c(-1, 1), na.rm = TRUE, theme.bw = TRUE)